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Abstract 

Background: Many Poaceae species show a gametopliytic selfMncompatibility (GSI) system, which is controlled by 
at least two independent and multiallelic loci, S and Z. Until currently, the gene products for S and Z were unknown. 
Grass SI plant stigmas discriminate between pollen grains that land on its surface and support compatible pollen tube 
growth and penetration into the stigma, whereas recognizing incompatible pollen and thus inhibiting pollination 
behaviors. Leymus chinensis (Trin.) Tzvel. (sheepgrass) is a Poaceae SI species. A comprehensive analysis of sheepgrass 
stigma transcriptome may provide valuable information for understanding the mechanism of pollen-stigma 
interactions and grass SI. 

Results: The transcript abundance profiles of mature stigmas, mature ovaries and leaves were examined using 
high-throughput next generation sequencing technology. A comparative transcriptomic analysis of these tissues 
identified 1,025 specifically or preferentially expressed genes in sheepgrass stigmas. These genes contained a significant 
proportion of genes predicted to function in cell-cell communication and signal transduction. We identified 1 1 1 
putative transcription factors (TPs) genes and the most abundant groups were MYB, C2H2, C3H, PARI, MADS. 
Comparative analysis of the sheepgrass, rice and Arabidopsis stigma-specific or preferential datasets showed broad 
similarities and some differences in the proportion of genes in the Gene Ontology (GO) functional categories. Potential 
SI candidate genes identified in other grasses were also detected in the sheepgrass stigma-specific or preferential 
dataset. Quantitative real-time PGR experiments validated the expression pattern of stigma preferential genes including 
homologous grass SI candidate genes. 

Conclusions: This study represents the first large-scale investigation of gene expression in the stigmas of an SI grass 
species. We uncovered many notable genes that are potentially involved in pollen-stigma interactions and SI mechanisms, 
including genes encoding receptor-like protein kinases (RLK), GBL (calcineurin B-like proteins) interacting protein kinases, 
calcium-dependent protein kinase, expansins, pectinesterase, peroxidases and various transcription factors. The availability 
of a pool of stigma-specific or preferential genes for L. cliinensis offers an opportunity to elucidate the mechanisms of SI 
in Poaceae. 
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Background 

Self-incompatibility (SI) is a widely distributed mechanism 
to prevent self-fertilization in flowering plants [1,2]. Clas- 
sic genetic studies revealed two major classes of homo- 
morphic SI, gametophytic (GSI) and sporophytic (SSI) 
systems [3]. Recent studies have focused on identifying 
molecules that mediate SI events. In Solanaceae, Rosaceae 
and Plantaginaceae GSI species, the stylar inhibition of 
incompatible pollen is mediated through an interaction 
between a stylar S-RNase (female determinant) and the 
pollen tube-borne F-box protein, SLF [4,5]. GSI in poppies 
(Papaveraceae) is mediated by a calcium-based signaling 
pathway that is most likely triggered by an interaction 
between a stigma-expressed S -glycoprotein ligand (a small 
extracellular signaling molecule) and its cognate plasma 
membrane receptor, which is a pollen-expressed Ca ^- 
channeled protein [6]. SSI in Brassicaceae is regulated 
by an S-haplotype-specific protein interaction between 
a stigma-specific S-receptor kinase (SRK, female deter- 
minant) and its cognate ligand S-locus cysteine-rich 
(SCR, male determinant) [7,8]. 

Poaceae exhibit a GSI system controlled by at least 
two independent and multiallelic loci (S and Z) [9,10]. 
An incompatibility reaction occurs when both the S and 
Z alleles of the haploid pollen are matched in the pistil. 
The mechanisms of the grass GSI are poorly understood 
relative to the well-characterized single-locus GSI system 
mentioned above. Numerous efforts to identify the genes 
of the S and Z loci have been performed, but currently, 
the genes for S or Z on either the pollen or pistil side 
remain unknown [10,11]. 

Li et al. [12] reported the putative pollen S-gene BM2 
identified from Phalaris coemlescem, which is a grass SI 
species. This gene was originally observed to co-segregate 
with the S genotype. Later studies indicated that BM2 en- 
codes a functional thioredoxin protein, which is possibly 
involved in post-translational protein modification and is 
located at ca. 1 cM from the S-locus [9]. Hackauf and 
Wehling [13] identified a putative ubiquitin-specific prote- 
ase (UBP) gene with a pistil-specific expression pattern and 
co-segregation with the Z-locus, thus suggesting the poten- 
tial involvement of protein ubiquitination in the grass SI 
system. However, whether the UBP gene is a component of 
the Z-locus or a linked gene with suppression recombin- 
ation around the Z locus is unknown [14]. Kakeda [15] iso- 
lated two F-box genes in Hordeum bulbosum based on the 
partial clone HAS175, which is an anther-derived marker 
showing complete linkage to the S locus. Yang et al. [11] 
constructed subtractive hybridization cDNA libraries to 
enrich the SI expressed genes in Lolium perenne. 

There are certain physiological characteristics of grass SI 
that resemble features of the Papaver system and Brassica 
SSI including the dry stigma type and a rapid response to 
stigma-pollen interactions [14]. In these three SI systems. 



the stigmas critically discriminate between compatible 
and incompatible pollen grains and promote compat- 
ible pollen growth while preventing incompatible pollen 
adhesion, hydration, germination or invasion. To fulfill 
these unique functions, stigmas should express genes that 
encode the proteins required for appropriate development 
of and accomplishment in pollen-stigma interactions. It 
is likely that many of these genes are stigma-specific or 
preferential. To gain a better understanding of the mo- 
lecular basis of stigma development and pollen-stigma 
interactions, it is important to generate a near-complete 
catalog of genes expressed specifically or preferentially 
in stigma cells. 

Currently, genome-wide transcriptome analyses of wet, 
dry and semidry stigmas have been conducted on tobacco 
[16], Arabidopsis [17,18], rice [19], maize [20] and Senecio 
squalidus [21]. However, stigma transcriptome data sets 
are not available in SI species of the Poaceae, the fourth 
largest family of flowering plants, which includes all 
cereal crops and major cultivated forage crops [22]. To 
develop insight into the conserved and diverse aspects 
of stigma development and functions across multiple 
species, it is useful to expand the available transcrip- 
tomic resources for incompatible species to allow com- 
parative analysis. 

Leymus chinensis (Trin.) Tzvel. (sheepgrass), a member 
of the tribe Triticeae, is a key dominant species of typical 
grassland communities in the Eurasian steppe region 
with economic and ecological value. The species is one 
of most important forage crops for grazing animals in 
China and is important in environmental protection, in- 
cluding soil and water conservation [23]. L. chinensis is 
a GSI species with incompatibility responses occurring 
on the stigma surface before or immediately after pollen 
germination [24]. 

The aims of the present study were to search for genes 
expressed specifically or preferentially in the stigma of the 
incompatible grass plant at the entire genome level using 
high-throughput next generation sequencing technology to 
perform an RNA-seq analysis on the stigmas, leaves and 
ovaries of L. chinensis. Here, the gene expression profiling 
of various tissues identified 1,025 genes that were predicted 
to be specifically or preferentially expressed in the stigma. 
The stigma gene sets were significantly enriched in cell- 
cell communication and signal transduction mechanisms, 
intracellular trafficking, secretion and vesicular transport. 
Our results also indicated that the sheepgrass stigmas 
shared common SI candidate genes, such as kinase genes, 
identified in Lolium. To our knowledge, this study is the 
first large-scale survey of gene expression in the stigmas 
of an SI grass species. The availability of a pool of stigma- 
specific or preferential genes for L. chinensis offers an 
opportunity to elucidate the molecular mechanisms of 
pollen-stigma interactions and SI in Poaceae. 
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Results 

The identification of stigma-specific or preferential genes 

To identify the genes involved in stigma development 
and required for unique stigma functions, we performed 
an RNA-seq analysis on the stigmas, leaves and ovaries 
of Leymus chinensis. These three samples enabled us to 
distinguish stigma-specific transcripts from transcripts 
that contribute to common cellular activities and functions. 
Mature stigmas, which were unpoUinated and reached the 
largest expansion, were collected from the pistils just be- 
fore anther dehiscence to avoid contamination with pollen 
grains (Figure 1). The stigmas were excised by cutting the 
pistil just below the base of the plumose stigma, and the 
remainder of the pistil was used as an ovary sample. Ap- 
proximately 1,500 stigmas, 1,500 ovaries and 20 leaves 
of 10-day old seedlings were collected for each RNA 
preparation. The mRNA from the three tissues was used 
to construct cDNA libraries, which were then sequenced 
on an lUumina HiSeq 2000 system. Paired-end 100-bp se- 
quence reads were then generated. After quality control, 
8,351,351, 6,046,077 and 5,654,202 high-quality reads were 
obtained from the stigmas, ovaries and leaves libraries, 
respectively. 

In previous studies, we performed long-read transcrip- 
tome sequencing of several cDNA libraries from differ- 
ent tissues of sheepgrass, including flowers before and 
after pollination and leaves, using a Roche 454 GS FLX 
sequencer [23]. A total of 87,214 unigenes were obtained, 




Figure 1 The mature pistil that was collected just before anther 
dehiscence. A long red line showing the position where a stigma 
was excised and a short red line showing the cutting position 
between ovary and rachis. 



including 32,416 contigs and 54,798 singletons. The mean 
contig size was 607 bp. To identify the genes that cor- 
responded to the sequencing reads in each library of 
the present study, the purity-filtered high-quality reads 
were mapped against the reference transcriptome data- 
set, which was generated using Roche 454 pyrosequencing 
technology. Altogether, 26,325, 29,748 and 30,620 contigs 
were detected at the stigmas, ovaries and leaves, respect- 
ively (Additional file 1). The stigma contigs, ovary contigs 
and leaf contigs represent about 81.2%, 91.8% and 94.4% 
of the total contigs of the Roche 454 sequencing dataset, 
respectively. Of these unigenes, 25,501 were common to 
all three tissues (Figure 2). A substantial number of genes 
showed tissue-specific expression: 183 in the stigma, 683 
in the ovary, and 1,432 in the leaves. 

Using a fold change of >1, false discovery rate (FDR) 
of <le-05 as a cutoff to identify differentially expressed 
genes, 842 contigs were up-regulated in the stigma 
compared to the ovary and leaf (Additional file 2 and 
Additional file 3). These contigs were considered to be 
preferentially expressed in the sheepgrass stigma. Thus, 
1,025 (183 plus 842) contigs were specifically or preferen- 
tially expressed in the stigma (hereafter designated the 
stigma dataset) (Additional file 4). 

Confirmation of stigma-specific or preferential genes by 
real time PCR analysis 

To validate the expression of the genes detected by the 
RNA-seq analysis, 18 genes with different expression 
patterns in the stigmas and leaves were selected for real- 
time PCR analysis (the primer sequences are available in 
Additional file 5). Scatterplots were generated by compar- 
ing the log2 fold change determined by the transcriptome 

( N 

Stigma (26325) 




Ovary (29749 ) Leaf (30620 ) 

Figure 2 A Venn diagram showing relationships between the 
three transcriptome datasets (stigmas, ovaries, and leaves). 

Numbers in parentheses are total number of expressed genes in 
each tissue type. 
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analysis and real time PCR. The correlation between these 
two analyses was then evaluated. The results showed that 
the expression patterns of these genes examined by real 
time PCR were well correlated with those by RNA-seq 
(R^ = 0.811), thus supporting the reliability of our RNA- 
seq data (Figure 3). 

KOG annotation of the stigma-specific or 
preferential genes 

From the 1,025 contigs specifically or preferentially 
expressed in the stigma, 648 contigs had putative func- 
tions and 377 contigs were assigned as hypothetical 
proteins (Additional file 4) based on BLASTX similarity 
searches against the Swiss-Prot and NCBI non-redundant 
protein sequences (nr) databases. The 648 contigs were 
further annotated and classified based on EuKaryotic 
Orthologous Groups (KOG) categories. Altogether, 522 
contigs were assigned KOG functional annotations 
(Additional file 4) and grouped into 23 functional cat- 
egories (Figure 4). Among the KOG categories, "Signal 
transduction mechanisms" (14.8%) was the most abun- 
dant (Additional file 6). 

Gene ontology (GO) annotation 

Gene ontology (GO) assignment programs were utilized 
for the functional categorization of annotated unigenes 
of the sheepgrass stigma dataset. In many cases, multiple 
terms were assigned to the shared contig. A total of 420 
contigs were categorized into 40 main functional groups 
belonging to 3 categories (biological processes, molecular 
functions, and cellular components). The GO terms of 
their subcategories are presented in Figure 5. 




♦ ♦-4 - 



-6 L 

Figure 3 Validation of RNA-seq results by quantitative real time 
RT-PCR. Correlation plots indicating the relationship between qPCR 
results (fold change; Y-axis) of 18 selected genes expressed in the 
stigmas and leaves and the corresponding data from RNA-seq 
analysis (X-axis). 



A comparison of the GO classification of stigma-specific 
or preferential genes in sheepgrass, rice and arabidopsis 

The stigma-specific or preferential expressed gene data- 
sets of rice [19] and Arabidopsis [17], two compatible spe- 
cies, are available to provide the opportunity to compare 
the diversity and conservation of stigma-specific or prefer- 
ential genes in species of mono- and dicotyledons. A total 
of 462 unigenes of rice stigmas and 115 of Arabidopsis 
stigmas were categorized into functional groups using the 
GO assignments. The WEGO online tool [25] was used to 
classify the GO functions for all genes of the sheepgrass, 
rice and Arabidopsis stigma-specific or preferential data- 
sets to understand the distribution of gene functions at the 
macro level. The comparative analysis showed that these 
three stigma transcriptomes shared broad similarities in 
the proportion of genes in the three main categories and 
many subcategories (Figure 5). However, differences were 
also detected among the three species. The high pro- 
portions of genes involved in the GO terms "protein 
tag", "translation regulators", "cell killing", and "rhyth- 
mic process" were only represented in the sheepgrass 
stigma dataset. The genes of the GO term "symplast" 
were only observed in Arabidopsis. These results suggest 
that despite the differences in the relationship and gross 
stigma morphology among these three diverse species, 
certain conserved mechanisms of structural development 
and the pollen-stigma interaction process appear to be 
common. Conserved stigma-specific genes may be evi- 
dence of the convergent evolution of classes of genes to 
acquire stigma-specific functions in diverse angiosperm 
species. Unique genes in the sheepgrass stigma dataset 
may imply intrinsic differences in stigma development 
and the pollination process between compatible and in- 
compatible species. 

GO enrichment analysis 

A singular enrichment analysis (SEA) [26] was performed 
to identify the significantly enriched GO terms in genes 
specifically or preferentially expressed in sheepgrass stigma. 
The results showed that 19 GO terms were overrepre- 
sented in the stigma based on a P-value < 0.001 and the 
FDR < 0.05 cutoffs, which were 13, 3 and 3 for the cellular 
components, biological processes, molecular functions 
categories, respectively. The detailed results of the Go 
enrichment analysis are presented in Table 1. 

Kyoto encyclopedia of genes and genomes (KEGG) 
pathway mapping 

The genes of the sheepgrass stigma dataset were anno- 
tated based on the Kyoto Encyclopedia of Genes and 
Genomes (KEGG) databases [27]. The results showed 
that KEGG numbers could be assigned to 128 contigs 
(Additional file 4). These unigenes were further annotated 
to BRITE functional hierarchies. The BRITE functional 
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Figure 4 KOG function classification. The unigenes were aligned to the KOG database to predict and categorize possible functions. A total of 
522 unigenes were assigned to 23 categories. 



mapping categorized the unigenes into diverse KEGG 
pathways (Figure 6). 

Putative transcription factors in the sheepgrass 
stigma dataset 

Understanding transcription factors (TFs) patterns of ex- 
pression is of particular interest because the expression and 
activity of these regulatory proteins is crucial for controlling 
the expression of numerous genes, and thus, their ability to 
regulate biological pathways and developmental processes. 
We searched the TFs in the sheepgrass stigma dataset using 
a BLASTX search against the PlnTFDB (Plant Transcrip- 
tion Factor Database) (version 3.0) (E-value < le-10). We 
identified 111 putative TF unigenes from the sheepgrass 
stigma dataset. These unigenes were clustered into 26 TF 
families. The top 19 TF families in the sheepgrass stigma 
dataset are shown in Figure 7. 

Sheepgrass stigmas share certain potential SI candidate 
genes identified in other grasses 

Yang et al. [11] identified 10 expressed SI candidate genes 
for the S and Z loci in L. perenne from the SI cDNA 



libraries. We compared gene sequences from Lolium with 
the sheepgrass unigene set to reveal the conservation of 
stigma preferential genes among different SI grass stigmas. 
The nucleotide sequences of 10 Lolium SI candidate genes 
(Canl35, accession number AM991118; CanlSO, AM99 
1119; CanlO, AM991120; Canl36, AM991121; Canl8, 
AM991122; Can3, AM991123; Can94, AM991124; Can4, 
AM991125; Canl39, AM991126; and Canl51, AM991 
127) were extracted from GenBank and blasted against 
the sheepgrass dataset containing 87,214 unigenes. Signifi- 
cantly similar sequences (using an E value < le-05) were 
identified for all 10 Lolium SI candidate genes. Nine of 
the ten best hits were included in the sheepgrass stigma 
dataset and four (contig23066, contigl5874, contig23698, 
and contig24204) were included in the 1,025 specific or 
preferential unigenes dataset (Table 2). An RT-PCR 
analysis confirmed the stigma-preferential expression 
for these four contigs (Figure 8). The Lolium SI candi- 
date genes that had the best hits in the sheepgrass 
stigma gene dataset included all three predicted kinase 
genes and a kinase-partner gene of the 10 SI candidate 
genes, i.e., Can3 (predicted as a CBL-interacting serine/ 
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Figure 5 GO assignment of all identified genes expressed specifically or preferentially in the stigmas of sheepgrass, rice and 
Arabidopsis. The genes were mapped to three main categories: biological processes, cellular components and molecular functions. The right 
hand y-axis indicates the number of annotated genes. 



threonine-protein kinase 15), Can94 (predicted as a 
calcium-dependent protein kinase), Can4 (predicted as 
a serine/threonine-protein kinase NAK), and Canl36 
(predicted as a pollen-specific kinase partner protein). 
Two (Can3 and Can4) of the three predicted kinase 
genes and the kinase-partner gene showed significant 
similarities to the unigenes of the specific and prefer- 
ential sheepgrass stigma gene set, thus indicating a po- 
tential role for kinase activity and conserved functional 
mechanisms in the grass SI response. 

It was proposed that Bm2 identified from P. coeru- 
lescens was involved in the SI response through the 
post-translational modification of other proteins by thior- 
edoxin proteins [9] . To detect whether the BM2 gene was 
expressed in sheepgrass, we blasted the gene sequence 
of BM2 (GenBank: AF159388.1) against the sheep- 
grass unigene set. The result indicated that Bm2 matched 
three unigenes (contigl2829, contigl2830, and contigl 
2832) with significantly similar sequences (identity > 80%, 
E-value < le-10). The RT-PCR experiments validated that 
contigl2830 was expressed in mature sheepgrass stigmas, 
thus suggesting that the BM2 gene plays an important role 
in the SI process. 

Discussion 

The present studies have identified 1,025 genes pre- 
dicted to be specifically or preferentially expressed in the 



stigmas of sheepgrass using high-throughput next gener- 
ation sequencing technology to perform an RNA-seq 
analysis. This is a minimium estimate because of strin- 
gent filtering process we used in analyzing HiSeq 2000 
sequencing data, incomplete detection of the expressed 
genes based on Roche 454 sequencing, decrease in the 
number of identifed genes mapped against the Roche 
454 reference dataset compared to those of actually ex- 
pressed genes in each tissue. It is likely that at least some 
of the specific or preferential sheepgrass stigma genes iden- 
tified in this study function in development of the stigma 
papillary cells and pollen-stigma interactions. Our analysis 
of the stigma transcriptome demonstrates the unique 
features of the stigma transcriptional profile. The func- 
tional specialization of the stigma for extracellular 
interactions is reflected by an increased number of 
stigma-specific or preferential expressed genes involved 
in signal transduction, such as kinases, receptor-like ki- 
nases and calcineurin B-like proteins (CBL), the manu- 
facturing of enzymes that most likely facilitate the entry 
of the pollen tube into the stigma, transcription factors 
that are possibly implicated in the pollination process, 
and defense-related genes that function in the defense 
against pathogens or potentially in response to pollin- 
ation. It is hypothesized that some of these genes are 
potentially involved in the pollination response and SI 
mechanisms. 
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Table 1 The overrepresented functional GO terms of the sheepgrass stigma-specific or preferential genes 


GO term 


Ontology 


Description 
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Cell-cell communication and signal transduction 

Genes predicted to function in cell-cell communication 
and signal transduction are particularly notable in the 
context of pollen-stigma interactions and SI. The sheep- 
grass stigma dataset contained a significant proportion 
of genes consistent with the primary recognition role of 
the stigma in discriminating genotypes of pollen grains. 
The most noteworthy are genes encoding various puta- 
tive protein kinases, including probable receptor-like 
protein kinases (RLK), mitogen-activated protein kinases 
(MAPKs), serine/threonine-protein kinase, CBL-interacting 
protein kinases, and calcium-dependent protein kinase. 

A stigma-specific receptor kinase, SRK, which was most 
extensively studied in Brassicaceae, interacts with the lig- 
and (SCR) held on the poUen coat and initiates a signal 
transduction cascade that inhibits the growth of self-pollen 
and the pollen tube [28]. In the Papaver GSI mechanism, 
the female S-determinant (PrsS) is a small stigma-secreted 
signal peptideand acts as a signaling ligand. The male 
S-determinant (PrpS) has a small predicted extracellu- 
lar domain that has been shown to be involved in both 
binding to PrsS and mediating SI [29]. A model for the 
operation of SI in Papaver has proposed that PrsS acts as a 
signaling ligand that interacts with PrpS. The sheepgrass 



stigma dataset contained 3 contigs, i.e., contig24113, 
contig28482, and contig24154, which are predicted to 
encode probable receptor-like protein kinase At2g42960, a 
probable LRR receptor-like serine/threonine-protein kinase 
Atlg63430, and somatic embryogenesis receptor kinase 1, 
respectively. A homologous gene encoding receptor-like 
protein kinase Atlg63430 was also identified in the Senecio 
squalidus SSI stigma preferential gene dataset. The con- 
tig24113, contig28482 and contig24154-encoded proteins 
may have a potentially important role in cell-cell commu- 
nication and signal transduction in the context of poUen- 
stigma interactions and are therefore notable candidates 
for further studies of sheepgrass stigma SI mechanisms. 

MAPKs are key contributors in signal transduction 
mechanisms and regulate important cellular processes, in- 
cluding cell proliferation, survival, and death in eukaryotes 
[30]. In addition, MAPKs are known to be involved in 
the activation of plant defense responses and result in 
programmed cell death (PCD) and pathogen resistance 
[30,31]. A MAPK, p56, is responsible for mediating SI 
induced PCD in Papaver [32,33]. p56 is involved in the 
loss of pollen viability, stimulation of caspase 3 like 
(DEVDase) activity and later DNA fragmentation in in- 
compatible pollen [33]. We also identified two contigs. 
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Figure 6 Histogram presentation of the KEGG pathway categories. The Y-axis indicates the number of unigenes assigned to a specific 
pathway. The X-axis indicates the KEGG pathway. 



contigl3945 and contig34021, in the sheepgrass stigma 
dataset that are predicted to encode for mitogen- 
activated protein kinase 10 and mitogen-activated pro- 
tein kinase 1, respectively. The presence of MAPKs in 
the sheepgrass stigma gene dataset impUes that they are 
Ukely involved in the early steps of pollination and SI 
responses. 

Temporal and spatial changes in cytosolic Ca^* concen- 
trations generate specific Ca ^ signals, which are decoded 
by Ca^^ sensors [34,35]. As a class of calcium-sensing 




Figure 7 Top 19 transcription factor families of tlie Leymus 
chinensis stigma dataset. The Y-axis indicates the number of 
unigenes assigned to a specific TF family. The X-axis indicates the 
top 19 TF families. 



proteins, calcineurin B-like proteins (CBLs) specifically 
target a group of sucrose non-fermenting-related serine/ 
threonine kinases (SnRK3), named CBL-interacting pro- 
tein kinases (CIPK), to mediate the sensed calcium signal 
[36-40] . The CBL-CIPK system is involved in a wide range 
of signaling pathways, including abiotic stress responses to 
drought and salt, innate immunity, plant hormone re- 
sponses and channel regulation [41,42]. No CBLs or 
CIPKs have been implicated as signaling components in 
plant SI mechanisms. It is particularly notable that four 
contigs encoding CBL were observed in the analysis of the 
sheepgrass stigma dataset. These novel genes were not 
previously identified as pistil-specific or preferential in 
other species. It was previously mentioned that the Can3 
(AM991123) gene, which putatively interacts with CBL, 
was identified in L. perenne stigma [11]. These results 
provide new insights into the function of CBL during 
pollen-pistil interactions in plants and suggest that CBLs 
are potentially involved in SI mechanisms. 

In addition, the sheepgrass stigma dataset also contained 
many notable unigenes that encode calcium-related pro- 
teins, such as calmodulin-like proteins, probable calcium- 
binding proteins (CML49, CML13, CML18, CML16, CML 
22, CML21), and calcium-dependent protein kinases. The 
role of cytosolic free Ca ^ in the regulation of pollen tube 
growth has been well characterized [43,44]. It has been 
shown that calcium signaling is involved in the self- 
incompatibility response of Papaver and Bmssica [45,46] . 
In analogy to Papaver and Brassica SI, there is some 
evidence for the involvement of Ca^^-mediated signal- 
ing in grass SI. The application of Ca * antagonists to 
isolated stigmas resulted in the inhibition or delay of 
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Table 2 The unigenes of Leymus chinensis transcriptome dataset that best hit SI candidate genes in Lolium perenne 



Can - homology to Rice- annotation 


The best hit to the Leymus chinensis unigene dataset 


Identity (%) 


P-value 


Beta-expansin 2 precursor 


contig34758 (not included in stigma) 


89.81 


5.00E-46 


Pathogenesis-related protein PRBl-3 precursor 


4-GJVU7SP04IOQ46 (in stigma) 


93.52 


3.00E-38 


Expressed protein 


contig26721 (in stigma) 


91.98 


2.00E-134 


Pollen-specific kinase partner protein 


contig23066 (included in 1025 unigene) 


91.52 


2.00E-125 


ELMO domain-containing protein 2 


contigl5874 (included in 1025 unigene) 


90.45 


5.00E-106 


CBL-interacting serine/threonine-protein kinase 15 


contig23598 (included in 1025 unigene) 


91.22 


0 


Calcium-dependent protein kinase 


contig24124 (in stigma) 


92.53 


0 


Serine/threonine-protein kinase NAK 


contig24204 (included in 1025 unigene) 


90.13 


3.00E-177 


gtkl6 protein 


contig15929 (in stigma) 


95 


1 .OOE-47 


myb-like DNA-binding domain, SHAQKYF class 


familyprotein contig 05419 (in stigma) 


92.31 


1 .OOE-36 



the SI response on self-pollination in S. cereale and 
L. perenne [10,47]. It is tempting to hypothesize that 
stigma S and Z determinants act as signal molecules 
in grass by interacting with the "self pollen S and Z 
partners at the tube plasma membrane, thus inducing a 
Ca ^-mediated signal cascade that results in pollen tube 
inhibition. However, whether a Ca ^-mediated signaling 
cascade is involved in grass SI is unknown. 

Many recent studies have also implicated ROS/H2O2 
as signaling molecules involved in plant reproductive 
processes such as pollen tube growth [48-50] and pollen- 
stigma interactions [51,52]. Because peroxidases can gener- 
ate and consume H2O2 [53], they should be considered as 
potentially important components of signal transduction 
pathways. Senecio squalidus stigmas accumulate high 
amounts of ROS, particularly H2O2, in their epidermal cells 
(papillae). The first plant peroxidase gene, SSP (stigma- 
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6 
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4 - 
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III 



■contig23066 
■contigl5874 

contig23698 
□contig24204 



leaf 



stigma ovary 

Figure 8 Expression pattern of four unigenes that best hit SI 
candidate genes of Lolium perrene. The relative expression levels 
are represented in arbitrary units (A.U.) normalized to the expression 
level of the actin gene, used as a reference, in each sample. 



specific peroxidase), expressed specifically in stigma papil- 
lae was identified in S. squalidus [54]. In the sheepgrass 
stigma dataset, we detected two contigs that encode per- 
oxidase. The potential signaling role of these interesting 
genes deserves further investigation. 

Enzymes most likely facilitating entry of the pollen tube 
into the stigma 

In flowering plants with dry stigmas, the stigmatic sur- 
face is covered by a continuous layer of cuticle, which is 
a significant barrier to pollen tube entry [55]. For suc- 
cessful pollination and subsequent fertilization, this bar- 
rier must be breached by cutinases from the pollen grain 
surface and stigma pellicle [56]. 

Cutinase is an esterase that degrades cutin [57-59]. 
Active cutinases have been detected in the pollen of 
Bmssica [60] and Tropaeolum [61], in which they pre- 
dominantly localize to the intine region of the pollen 
wall. Knox et al. [62] has suggested that esterases from 
the pollen and stigmatic pellicle combine to form an 
"active cutinase complex" which degrades the cuticle in 
the region of penetration. Although cutinases have been 
implicated in the process of pollen tube penetration in 
many species with dry stigmas, none have been fully 
characterized at the molecular level [63]. The identifica- 
tion of the genes that encode plant cutinases is important 
for elucidating the digestion mechanism of the stigmatic 
cuticle during pollen-stigma interactions. The sheepgrass 
stigma dataset contains three contigs that were predicted 
to encode GDSL esterase, two contigs for probable pectin 
esterase, and four contigs for pectin esterase inhibitor. 
These unigenes are very interesting candidates for further 
studies on enzymatic degradation of the stigmatic cuticle 
during pollination. 

When the cuticle has been breached, pollen tubes in- 
vade the stigmatic papillary cell wall by growing between 
its outer and inner layers. The pollen tube then con- 
tinues to grow through the intercellular spaces (extracel- 
lular matrix, ECM) of the transmitting tissue towards 
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the ovary [64,65]. Much evidence has indicated that 
pollen- and stigma-derived enzymes and proteins most 
likely facilitate the entry of the pollen tube into the 
stigma by loosening the papillary cell wall and the ECM 
and/or modifying the growing pollen tube wall. These 
enzymes and proteins included expansins [66], pectin 
esterases [67], polygalacturonases [68], glucanases [69], 
xylanases [70], pectate lyases, pectin methylesterases [67], 
and peroxidases [71]. Lipid transfer proteins (LTP), one of 
the most abundant proteins present in the stigma of sev- 
eral species, have cell wall-loosening activity that is similar 
to expansion, despite its assumed role in lipid metabolism. 

The sheepgrass stigma dataset contained many unigenes 
for putative cell wall-localized proteins that might func- 
tion in regulating the expansion or loosening of the stig- 
matic papillary cell wall and ECM, thus regulating the 
entry of the pollen tube into the stigma. Among these pro- 
teins are two expansins, five pectinesterase, a pectin ester- 
ase inhibitor, two peroxidases, five polygalacturonases, 
three pectate lyases, and seven lipid transfer proteins. 
Because their better studied relatives exhibit cell wall- 
loosening activity and have been implicated in ECM 
remodeling and the disassembly of ECM structural com- 
ponents, some of these proteins (e.g., expansins and hydro- 
lase) may facilitate the penetration of the pollen tube into 
the stigma by modifying the papillary cell wall and ECM. 

Expansins are cell wall-loosening proteins that regu- 
late cell wall expansion and cell enlargement in a pH- 
dependent manner [72]. The expansin superfamily (EXP) 
comprises four distinct families: expansin A (EXPA), expan- 
sin B (EXPB), expansin-like A (EXLA) and expansin-like B 
(EXLB) [73]. Experimental evidence has indicated that 
EXPA and EXPB proteins function in the loosening and 
extension of cell walls, whereas the exact functions of 
EXLA and EXLB remain unknown [74,75]. The sheep- 
grass stigma dataset contained A and B expansins. The 
latter have been suggested to act as cell wall-loosening 
agents to facilitate pollen adhesion and pistil penetration 
[76]. Yang et al. [11] also identified a prominent SI can- 
didate, Canl35, representing a putative EXPB protein in 
L. perenne stigmas. Because the assumed functions of 
expansins have been implicated in the process of pollen 
tube penetration, it is possible that the expansins expressed 
preferentially in sheepgrass stigmas might be related to 
pollen-stigma interactions and SL 

Transcription factors possibly implicated in the 
pollination process 

The processes underlying pollen recognition and pollen 
tube growth require the concerted action of genes that 
are regulated by transcription factors [77]. A total of 111 
unigenes encoding the transcription factors of 26 fam- 
ilies were identified in mature sheepgrass stigmas. Most 
TPs (58%) belonged to seven famUies: MYB, C2H2, C3H, 



PARI, MADS, AP2-EREBP and HB. Among these fam- 
ilies, the MYB genes form the largest category. Most plant 
MYB genes encode proteins of the R2R3-MYB class [78]. 
Numerous R2R3-MYB proteins are involved in the con- 
trol of plant-specific processes including: primary and sec- 
ondary metabolism, cell fate and identity, developmental 
processes, and responses to biotic and abiotic stresses 
[79,80]. For example, several R2R3-MYB proteins encoded 
by AtMYBO/GLl, AtMYB23 and AtMYB66/WER are 
involved in the determination of epidermal cell type. 
AtMYBO and AtMYB23 control trichome initiation in 
shoots, whereas AtMYB66 controls root hair patterning 
[81]. Two closely related R2R3-MYBs, AtMYB88 and AtM 
YB124/FLP regulate stomatal differentiation and pattern- 
ing. AtMYB98 regulates synergid cell differentiation during 
female gametophyte development, pollen tube guidance 
and the formation of the filiform apparatus [82]. Canl51, 
which is an important SI candidate gene that represents a 
putative Myb-like protein, was identified in L. perenne stig- 
mas [11]. MYB transcription factors were also detected in 
the stigmas of many other species, such as Arabidopsis 
[17,18], rice [19], tobacco [16], maize [20], and crocus [83]. 
Similar to our studies, the transcriptome analysis of ma- 
ture maize stigmas also indicated that the MYB genes rep- 
resented the largest category of transcription factors. The 
preferential expression of MYB family genes in the mature 
or pollinated stigma of many species implicates its pos- 
sible role in the pollination process and pollen-stigma in- 
teractions. Recently, Gomez-Gomez et al. [84] isolated 
CsMYBl, which is a transcription factor belonging to the 
R2R3 family gene, from Crocus sativus. CsMYBl is highly 
expressed in stigma tissues and poorly expressed in tepals, 
whereas no transcript was detected in either anthers or 
leaves. CsMYBl's expression is developmentally regulated 
with no transcript detected in early stage stigmas but high 
levels occurring in later stages. It was assumed that 
CsMYBl is possibly involved in the control of stigma mor- 
phological development. Although MYB genes have been 
identified and isolated in the stigma of many species, the 
precise functions of these ubiquitous stigma TPs in the 
pollination process are uncharacterized. 

Defense-related genes 

The stigmas of flowering plants are stubbornly resistant 
to pathogen invasion. Both wet and dry stigmas possess an 
effective pathogen defense system [1] and express a wide 
range of pathogenesis-related genes that presumably con- 
tribute to this stigmatic defense system [17,18,85]. These 
defense-related genes might function not only in defense 
against pathogens, but also in response to pollination. The 
mechanisms originally evolved in defense against pathogens 
may have been recruited to identify and reject self-pollen in 
an SI response [1]. The evidence indicated many striking 
similarities in the process of epidermal cell penetration by 
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fungal spores and the pollen tube, including parallels in cell 
wall modification, the presence of phenolic derivatives, the 
accumulation of high levels of extracellular calcium, and 
the synthesis of pi-3 glucan callose [1,19]. 

Our transcriptome analyses have identified several 
stigma-specific or preferential genes that encode pro- 
teins potentially involved in defense responses, such as 
dihydroflavonol-4-reductase (DFR) and Bax inhibitor 1 
(BI-1). DFR are key regulatory enzymes of flavonoid bio- 
synthesis [86]. Flavonoids are ubiquitous plant secondary 
metabolites and are accumulated in different organs and 
tissues of plants [87]. Flavonoids are mainly involved in 
protecting plants against predation and pathogens (bac- 
teria and fungi) [88] . High amounts of stigmatic flavonoids 
may be toxic to potential pathogens but are not toxic to 
pollen. It can be speculated that an unknown signaling 
mechanism identifies a pollen grain as a pollen grain and 
not a fungal spore or bacteria. 

Plant cells can respond to various stimuli, including 
fungal toxins and biotic and abiotic stresses, by initiating 
programmed cell death. Bax inhibitor 1 proteins act as 
negative regulators of stress-related programmed cell death 
in plants and animals [89]. The stable overexpression of 
this cell death inhibitor in barley reduces susceptibility to 
the necrotrophic fungus Fusarium graminearum [90]. The 
transient or stable overexpression of BI-1 strongly supports 
the penetration of Blumeria graminis f. sp. hordei into bar- 
ley epidermal cells [90,91]. Here, we present the first report 
that genes encoding BI-1 proteins were expressed in the 
stigma of flowering plants. However, further analysis is 
necessary to better understand the involvement of these 
notable proteins in pathogen defense and pollen-stigma 
interactions. 

Conclusion 

For the first time, we present a large-scale investigation of 
gene expression in the stigmas of sheepgrass (L. chinensis), 
a Poaceae GSI species, using a high-throughput RNA-seq 
analysis. Altogether, 26,325 contigs were detected in the 
stigmas. Of these contigs, 1,025 showed specific or prefer- 
ential expression. Quantitative real-time PGR confirmed 
the expression patterns of the genes examined by 
RNA-seq. This large-scale survey represents an effi- 
cient approach to identify the genes involved in the stigma 
pollination process. The functional specialization of the 
stigma for pollination behavior is reflected by unique fea- 
tures of its transcriptional profile. Many identified stigma- 
specific or preferential genes were potentially involved in 
encoding signaling molecules, such as protein kinases and 
receptor-like kinases, proteins and enzymes modifying the 
cell wall, cuticle or ECM to facilitate the entry of the 
pollen tube into the stigma, and defense-related proteins 
that function to defend against pathogens or in response 
to pollination. Most of the genes identified by the RNA- 



seq analysis were not previously studied or reported to be 
expressed in stigmas, and their potential involvement in 
stigma development and poUen-stigma interactions opens 
new avenues to better understand the grass GSI system in 
terms of function and evolution. The future functional 
characterization of these genes promises to elucidate the 
mechanisms that underlie incompatible pollen-stigma in- 
teractions in grass GSI systems. 

Methods 

Plant materials, RNA extraction, library construction and 
illumina sequencing 

The plant materials were harvested from sheepgrass plants 
grown under natural conditions in a field of the Beijing 
Botanical Garden, Chinese Academy of Sciences, Beijing, 
China. The stigmas were collected 1-10 hours before 
floret flowering by cutting the pistil below the base of the 
stigma, and the remainder of the pistil was used for the 
ovary samples. The leaves were collected from 10-day-old 
seedlings. 

The total RNA from the mature stigmas, mature ovar- 
ies and leaves was extracted using TRIzol Reagent ac- 
cording to the manufacturer's instructions (Invitrogen, 
Carlsbad, CA, USA). The quality of the total RNA was 
determined using a NanoDrop 2000 (Thermo Fisher, 
USA). The mRNA was purified from the total RNA sam- 
ples using a Dynabead mRNA Purification Kit according 
to the manufacturer's instructions (Invitrogen, Carlsbad, 
CA, USA), and the quality was assessed using an Aligent 
2100 Bioanalyzer (Agilent Technologies, Inc., Waldbronn, 
Germany). Double-stranded cDNA was synthesized using 
the Superscript Double-Stranded cDNA Synthesis Kit 
(Invitrogen, Carlsbad, CA, USA). Specific adapters were 
ligated to the fragmented cDNA and denatured to gen- 
erate single-stranded cDNA followed by emulsion PGR 
amplification. The sequencing was performed using an 
Illumina HiSeq 2000 sequence analyzer at Hanyu Genomics 
Institute (Shanghai, China). 

The sequence reads generated from each tissue sample 
were aligned to the reference transcriptome dataset ob- 
tained by Roche 454 pyrosequencing technology using 
SOAP2 software [92]. For the specific mapped uni- 
genes with uniquely matched reads, the transcript 
abundance level was normalized using the RPKM 
(Reads Per kb per Million reads) method. The RPKM 
values were computed as proposed by Mortazavi et 
al. [93]. If there was more than one transcript for a 
gene, the longest one was used to calculate its tran- 
script abundance and coverage. 

Real-time PCR analysis 

The total RNA (4 |ig) was reverse-transcribed with an 
oligo (dT) primer for cDNA synthesis using a Superscript 
III First-Strand Synthesis Kit (Invitrogen). Gene-specific 
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Additional file 3: Comparison of gene expression between the 
stigmas and ovaries of Leymus chinensis. 

Additional file 4: 1,025 genes specifically or preferentially 
expressed in the stigma of Leymus chinensis. 

Additional file 5: A list of real- time PCR primers used in the 
present study. 

Additional file 6: The Leymus chinenesis stigma specifically or 
preferentially expressed genes predicted to function in signal 
transduction mechanisms. 



primers were designed using PRIMEREXPRESS software 
(Applied Biosystems). The primer sequences are listed in 
Additional file 5. Quantitative PCR assays were performed 
in triplicate using SYBR Green Real-time PCR Master Mix 
(Toyobo, Osaka, Japan) with a Bio-Rad CFX96 Real-Time 
Detection System. The quantitative variation in the differ- 
ent replicates was calculated using the delta-delta thresh- 
old cycle relative quantification method. Amplification of 
L. chinensis actin was used as an internal control to 
normalize all data. 

Annotation, functional classification and pathway analysis 

All mapped contigs were annotated with GetORF from 
the EMBOSS package [94]. The ORE of each predicted 
protein was used for BLASTP searches against the Swiss- 
Prot and NCBI nr databases with thresholds of E-value < 
le-5. Domain-based alignments were performed against 
the KOG database at NCBI with a cut-off E-value of <le-5. 
GO annotations for describing the biological processes, 
molecular functions, and cellular components were ana- 
lyzed by GoPipe using a BLASTP search against the Swiss- 
Prot and TrEMBL databases with an E-value < le-5 [95]. 
The WEGO online tool [25] was used to classify the GO 
functions for all genes of sheepgrass, rice and Arabidopsis 
stigma-specific/preferential datasets to understand the 
distribution of gene functions at the macro level. An- 
other GO analysis tool, SEA [26], was used to identify 
overrepresented GO terms. KEGG pathways annotations 
were performed using the KEGG Automatic Annotation 
Server (KAAS) with the bi-directional best-hit information 
method [27]. KAAS annotates every submitted sequence 
with KEGG orthology (KO) identifiers that represent an 
orthologous group of genes directiy linked to an object in 
the KEGG pathways and BRITE functional hierarchy 
[27,96]. To identify transcription factor families represented 
in the sheepgrass stigma dataset, unigene sequences were 
searched against a complete list of transcription factor 
protein sequences of the Plant Transcription Factor Data- 
base (PlnTFDB: http://plntfdb.bio.uni-potsdam.de/v3.0/ 
downloads.php) using BLASTX with an E-value cutoff 
of < le-10. 

Availability of supporting data 

Sequence data from this article were deposited in the 
NCBI Sequence Read Archive (SRA065691). 

Additional files 



Additional file 1: Genes expressed in the stigmas, leaves, and 
ovaries of Leymus chinensis. RPKM: Reads Per Kilo bases per Million 
reads; NA: Not applicable. 

Additional file 2: Comparison of gene expression between the 
stigmas and leaves of Leymus chinensis. 
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